Insulator to superfluid transition in coupled photonic cavities in two dimensions 
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A two-dimensional square lattice of coupled photonic cavities is systematically investigated using 
quantum Monte Carlo simulations. The ground state phase diagram contains insulating phases with 
integer polariton densities surrounded by a superfluid phase. Finite-size scaling of the superfluid 
stiffness is used to extract the phase boundaries. The critical behavior is that of the generic, density- 
driven insulator-superfluid transition with dynamic exponent z — 2, with no special multicritical 
points with z — 1 at the tips of the insulating-phase lobes (in contrast to the Bose-Hubbard model). 
This demonstrates a limitation of the description of polaritons as structureless bosons. 
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One of the currently most active and intriguing areas of 
condensed matter physics is the search for "exotic" quan- 
tum phases and related quantum phase transitions |l[ in 
strongly correlated systems. However, the control of the 
relevant system parameters, which correspond to various 
coupling constants of theoretical model systems, remains 
restricted in chemically synthesized materials. This has 
stimulated investigations of artificially engineered nano- 
structures, in which a wider range of model behavior can 
be realized. Of these systems, coupled Josephson junc- 
tions and ultracold atoms in optical lattice are perhaps 
the most prominent so far 0, HI, 0j- The latter case is, 
among other things, ideally suited for implementing the 
Bose-Hubbard model H throughout its parameter space. 

Recently, another engineered structure was proposed 
which shows great promise for investigating novel quan- 
tum states of matter; arrays of optical cavities in each of 
which there is one or several atoms interacting with pho- 
tons according to the famous Jaynes-Cummings interac- 
tion [fj, |7|, l|| . The cavities are coupled through photon 
tunneling (hopping) between adjacent cavities. A par- 
ticular advantage of this structure over atoms in optical 
lattices is that the distance between the cavities is com- 
parably large which allows for the manipulation and 
measurement of individual cavities experimentally. In 
these systems, the number of photons is not conserved, 
but the number of polaritons — quasiparticles constructed 
as combinations of photons and atomic excitations — are 
conserved and have been argued to be analogous to the 
particles of the Bose-Hubbard model. Theoretically, it 
has been well established that there exists a polariton 
localization-delocalization transition 0, @, H, m, [Hi [H[ , 
in analogy with the insulator-superfluid transition in the 
Bose-Hubbard model Extending the analogy, a sys- 
tem with two species of photons has also been investi- 
gated and related to the two-component Bose-Hubbard 
model (l2j |. and a super-counter- fluidity phase has been 
predicted for circularly polarized photonic cavities [l3| . 
While experimental work on coupled optical cavities has 
not yet reached the stage where these states could be 



explored, Bose-Einstcin condensation of polaritons has 
been demonstrated in single two-dimensional cavities 
with a large number of embedded quantum wells . 

In this paper we will focus on the simplest version of 
optical cavities on the two-dimensional (2D) square lat- 
tice; in each cavity there is a single atom with two energy 
levels and the photons are polarized so that we need to 
consider only one-component photons. This system is 
described by the Hamiltonian 

H = -t^2(ala j + h.c.) + J2h i , (1) 

(ij) * 

where a, and a\ are photon annihilation and creation 
operators, and the single-site terms hi are given by 

hi = e(S? + ±) + uja\a % - /i(Sf + \ + at a*) 

+ g(S+a i + Sr a i). (2) 

Here a two-level atom is represented by a spin- 1/2, with 
being the corresponding raising and lowing opera- 
tors. The atomic excitation energy and the resonance 
frequency of the cavity are denoted by e and u>, respec- 
tively, and g is the coupling constant between photons 
and atoms. The excitations of isolated cavities (t = 0) 
are not photons but entangled photon-atom states — the 
polaritons. Their total number operator is given by 

Af = £< ol = 5>ja i + S? + ±), (3) 

i i 

which commutes with Hamiltonian, [J\f,7i] = 0, for all t. 
Thus fi in ([2]) is the chemical potential for polaritons. 

Exact numerical calculation with the density-matrix- 
renormalization-group (DMRG) method for the ID case 
has confirmed a ground state phase diagram with insu- 
lating and superfluid phases [lOl. A mean-field theory 
showed a similar behavior in 2D [7j. In both ID and 2D, 
the insulating phases have integer polariton density and 
form a series of lobes in the (t, fj,) plane, in qualitative 
agreement with the Bose-Hubbard model. Very recently, 
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Aichhorn et al. presented ID and 2D phase diagrams, as 
well as single-particle spectra, obtained with a systemati- 
cally controllable variational cluster approach ll[ . Their 
ID results are in good agreement with the DMRG calcu- 
lations. 

The previous numerical calculations have convincingly 
demonstrated the existence of the insulating and super- 
fluid phases and located their phase boundaries, but the 
exact nature of the phase transition has not yet been de- 
termined. A major difference between the cavity model 
(Q]) and the Bose-Hubbard model is that there is no direct 
on-site repulsion between photons or polaritons, which 
plays a key role [5] in the insulating phase of the Bose- 
Hubbard model (although there is an effective repulsion 
between photons mediated by the spins [15j|). Perhaps 
more importantly, the polaritons are complex quantum 
objects different from pure bosons. It is therefore not a 
priori clear that the phase transition between the insu- 
lating and superfluid phases should be in the same uni- 
versality class as in the Bose-Hubbard model. 

Here we investigate the quantum phase transition nu- 
merically, using approximate-free quantum Monte Carlo 
simulations (stochastic series expansion with loop up- 
dates [Hj]). We use finite-size scaling to extract the phase 
boundaries, as well as exponents governing the quantum- 
critical scaling behavior, and conclude that there is in- 
deed a qualitative difference between the two models. In 
the Bose-Hubbard model there are special multi-critical 
points at the tips of the insulating lobes, where the tran- 
sition occurs at constant density and the critical behav- 
ior is non-generic, with 3D classical XY exponents — most 
notably the dynamic exponent z = 1 at the tips. In con- 
trast, the transition in the optical cavity system is always 
of the generic, density-driven mean-field kind, with dy- 
namic exponent z = 2. Another interesting aspect of 
the phase diagram is the shape of the lobe tips, which 
are cusped in ID [lfj. The calculations of Aichhorn et 
al. also indicate cusps in 2D [ll| , although much less pro- 
nounced than in ID, in contrast to the smooth tips in the 
2D Bose-Hubbard model [B|, [13] • However, our calcula- 
tions show that the lobes are smooth. 

We consider two kinds of phase diagrams: (i) at chemi- 
cal potential p = in the (A,i) plane, where A = e — lo is 
the detuning parameter, and (ii) in theju, t) plane with 
the detuning A = (as studied also in [ll(). For simplic- 
ity we set g = e = 1. When the hopping t — 0, the Hamil- 
tonian can be readily diagonalized and we can extract the 
number of polaritons n po1 > as a function of the other 
parameters. For finite t, only the total number of polari- 
tons Af is conserved, and the average number (n po1 ) is in 
general non-integer. However, it is integer in the insu- 
lating phases at temperature T = 0. In our QMC simu- 
lations we have to impose a maximum polariton number 
n max P er site. The cut-off error is systematically checked 
by comparing results for n^l x = 3,4,5,6. We consider 
only the insulating phases with n po1 = 1,2 and find no 
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FIG. 1: (Color online) Superfluid density versus hopping for 
L = 14 (A), 18 (0), 20 (□), and 22 (O), at (/x, A) = (0,0.3). 
Error bars are at most comparable with the symbol size. The 
data extrapolated to L = oo, using an L~ 2 finite-size depen- 
dence, are shown (x) along with a fitted line. The inset shows 
the finite-size behavior corresponding to z = 2. 

detectable effects of the cut-off once nP°l K > 4. 

Using the finite-size-scaling hypothesis for the super- 
fluid density as a function of L, (3 = 1/T, and 5 = t — t c 
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we can determine a critical insulator-superfluid transi- 
tion point accurately by measuring p s (in the standard 
way through the photon winding number fluctuations) 
near the critical point for different lattice sizes L x L and 
hopping t. If we fix (3/L z = a, then the scaling function 
p() is a function of the single parameter SL 1 ^ . At the 
critical point 6 is zero and p s L d+z ~ 2 is then independent 
of lattice size, which implies that calculated p s L d+z ~ 2 
curves versus t for different L should intersect at the 
critical point. In our analysis we test z = 2, the generic 
dynamical exponent of Bose-Hubbard model, and, as al- 
ready mentioned, find this to hold on all points of the 
phase boundaries. We use two different aspect ratios, 
a — 1/3 and 1/2, for added confidence in the results. 

We first show in Fig. [T] examples of the behavior of the 
superfluid density close the phase transition. Results for 
the superfluid density in the range t = 0.0350 — 0.0382 
at fixed (/i, A) = (0.0, 0.3) are graphed for lattice sizes 
L = 14, 18, 20, and 22. Extrapolating the results to L = 
oo (using oc L~ 2 finite-size corrections, which our data 
exhibit), we observe the superfluid density in the thermo- 
dynamic limit becoming non-zero for t > t c = 0.03702(5). 
Above this critical point the superfluid density grows lin- 
early with 5, which if z = 2 is consistent with the gener- 
alized Josephson relation p s ~ S u i d + Z -' 2 ) in combination 
with the hyperscaling relation v = 1/z [5j. 

The finite-size scaling of the superfluid density at crit- 
icality, p s oc L~ z with z = 2, is shown in the inset of 
Fig. [2 The critical hopping t c is readily determined 
from the crossing point (more accurately than using the 
extrapolated p s above). We use this procedure to deter- 
mine the full phase boundaries of the first two insulat- 
ing lobes in the (i, A) plane at p = 0, shown in Fig. [5] 
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FIG. 2: (Color online) Phase boundary between the insulating 
and superfluid ground states in the (t, A) plane for polariton 
chemical potential fi = 0. The squares show two points deter- 
mined by finite-size scaling at fixed t, varying A. The other 
points were extracted by varying t at fixed A. The dot-dashed 
line shows the (n po1 ) = 1 polariton density contour. 



For comparison we also determine some critical points 
by scanning along fixed-t lines across the first insulating 
lobe (instead of scanning t at fixed A). As seen in the 
figure, the phase boundaries obtained by the two meth- 
ods are in excellent agreement. The lobe tips are lo- 
cated at (i, A) = (0.03872(5), 0.260(5)) for nP° l = 1 and 
(i, A) = (0.00358(4), 0.636(3)) for n po1 = 2. The t at the 
first lobe tip is about 11 times larger than at the second 
one, and, thus, the lobes shrink more quickly than the 
1/n behavior of the Bose-Hubbard model. This had al- 
ready been observed in the ID case 

[HI EH- Fig. [2] also 
shows the constant-density (n po1 ) = 1 contour inside the 
superfluid phase. Here this curve does not connect to 
the tip of the lobe, as it would have to do (as a require- 
ment for positive compressibility [B[) in a (t, //) phase 
diagram — we will discuss such a case below. In the Bose- 
Hubbard model, the universality class of the transition 
is different exactly at the tip of the lobe, because den- 
sity fluctuations are suppressed there due to an emergent 
particle-hole symmetry. This special transition therefore 
belongs to the three-dimensional XY universality class, 
whence 0=1. However, here we find z = 2 throughout, 
as we will comment on more extensively further below. 

Hyperscaling requires that the correlation length ex- 
ponent v = 1/z, and thus we should have v = 1/2. The 
onset superfluid density shown in Fig. [T] already demon- 
strated consistency with z = 2 and v = 1/2. We can test 
these exponents to higher precision by a standard data 
collapse procedure, graphing p s L 2 versus SL 1 /" for sev- 
eral system sizes. A typical example demonstrating the 
expected scaling is shown in Fig. [3) 

Aichhorn et al. considered a slightly different phase di- 
agram from the one we have focused on above. They 
calculated the phase boundary in the (t, p,) plane using 
a variational cluster method llj, which should become 
exact in the limit of large embedded clusters. In order to 
compare one-to-one with their results, we also extracted 
the boundary of the first insulating lobe in the same pa- 
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FIG. 3: (Color online) Finite-size data collapse of p„L z as a 
function of {t — t c )L}^ v , with z = 2 and v = 1/2, showing 
consistency with with the hyperscaling relation v = 1/z. 

rameter space, using the same finite-size scaling methods 
discussed above. Our results close to the first lobe tip 
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are shown in Fig. [4] along with the data of Ref. 
for different cluster sizes. At small hoppings, the phase 
boundaries obtained with the two methods agree very 
well, and there is very little dependence on the cluster 
size in the variational calculation. However, close to the 
lobe tip there are notable differences. Here the depen- 
dence on the cluster size in the variational calculation 
is much more pronounced, but our phase boundary does 
appear to be roughly consistent with the behavior for in- 
creasing cluster size. However, there is a qualitative dif- 
ference in that we obtain a smooth lobe tip, instead of the 
cusps that appear persistently for all cluster sizes in the 
variational calculation. It is quite possible that the cusps 
scale away with increasing cluster size, but carrying out 
this finite-size scaling would be difficult in practice, due 
to the limited range of available cluster sizes. We have 
checked our calculations at the lobe tip extensively, car- 
rying out the finite-sizse scalings using horizontal as well 
as vertical parameter scans, looking for possible drifts in 
the p s L 2 crossing points, and checking the dependence on 
the cut-off fim°ax- We do not find any indications of a cusp 
and conclude that the cusps are artifacts of small cluster 
sizes in the variational calculation. The figure also shows 
the (n po1 ) = 1 density contour in the superfluid phase, 
which connects to the tip of the lobe, as it should [B| . 

The inset of Fig. [4] shows the finite-size scaling of the 
superfluid density very close to the lobe tip. The location 
of the tip is determined to t = 0.05201(5), /j, = 0.185(5). 
Also here the scaling with z = 2 is excellent — there are no 
indications of drifts of the crossing points with increasing 
lattice size. Thus, we conclude that there is no special 
z = 1 multi-critical point. Although this different uni- 
versality, if present, should only exist at a single point, it 
should be reflected asaz = l^z = 2 finite-size cross- 
over in an extended region surrounding that point. In 
the Bose-Hubbard model z — 1 scaling close to the lobe 
tip can be seen ver y cl early for the same range of system 
sizes studied here [rjj (larger sizes were studied in this 
paper, but z = 1 scaling is seen already for L — 10). The 
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FIG. 4: (Color online) First-lobe boundary in the (t, fj,) plane 
for A = 0. The open circles show results of finite-size scaling 
using scans of t at fixed while the two solid squares were 
obtained at fixed t, scanning [i. The other curves are from 
the variational cluster calculations in Ref . p"H , with increa sing 
cluster size corresponding to the curves right to left (see [ill ] 
for details). The dot-dashed line touching the lobe tip is the 
(n po1 ) = 1 contour. The inset demonstrates scaling of the 
superfluid density with 2 = 2 close to the tip of the lobe, 
using data for L = 16, 18, 20, 22. 



tails of the localization mechanism are likely irrelevant 
(as long as the photon density fluctuations are not sup- 
pressed, which would lead to XY criticality) . There has 
been recent work based on a very small number (2-6) 
cavities to elucidate the nature of the effective particles 
in different parameter regimes of the coupled op- 

tical cavities. Two types of superfluids, photonic and 
polaritonic, were discussed 18( . However, it remains to 
be clarified whether there really is any clear distinction 
between the two, in particular as regards the phase tran- 
sition. In principle the method we have used here could 
also be used to probe the mechanism of localization on a 
more microscopic level, as a complement to studying the 
resulting macroscopic critical behavior. 
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to thank the ISSP for support and hospitality during an 
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Institute for Solid State Physics, University of Tokyo. 



fact that we observe essentially perfect scaling with z = 2 
everywhere is therefore very strong evidence in favor of a 
generic, mean-field density-driven transition also at the 
tip of the lobes. 

In summary, we have studied the quantum phase tran- 
sition separating a localized (insulating) and delocalized 
(superfluid) state in a 2D system of coupled optical cav- 
ities. Our main result is that this transition always has 
dynamic exponent z — 2, corresponding to the generic 
mean-field insulator-superfluid transition of interacting 
bosons. There are no special multi-critical points with 
z = 1, which do exist in the Bose-Hubbard model at the 
tips of the insulating-phase lobes. The multi-criticality, 
of the XY universality class, in the Bose-Hubbard model 
is a result of suppression of density fluctuations upon ap- 
proaching the lobe-tips [BJ . In the coupled optical cavities 
the polariton density fluctuations should also freeze out 
at the lobe tips, but the photon number is not conserved. 
The superfluid stiffness that we have studied here is de- 
fined soley in terms of the photons (calculated using their 
winding number fluctuations), and thus it appears that 
the superfluid is a state essentially describable in terms 
of the photons only (with the atomic excitations acting 
effectively as a fluctuating background potential). Al- 
though the effective particles in the insulating state are 
clearly polaritons, it is plausible that the phase transition 
can be described solely in terms of the photons, because 
although their number fluctuates in the insulating phase 
the essence of the transition is just that they become lo- 
calized. Our results that z = 2 everywhere supports this 
notion of a generic photon delocalized-localized transi- 
tion. On a microscopic level, the localization still of 
course involves the formation of polaritons, but the de- 
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